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We study the thermodynamic and structural properties of the superconducting vortex system in 
high temperature layered superconductors, with magnetic field normal to the layers, in the presence 
of a small concentration of strong random point pinning defects via numerical minimization of a 
model free energy functional in terms of the time-averaged local density of pancake vortices. Working 
at constant magnetic induction and point pinning center concentration, we find that the equilibrium 
phase at low temperature (T) and small pinning strength (s) is a topologically ordered Bragg glass. 
As T or s is increased, the Bragg glass undergoes a first order transition to a disordered phase which 
we characterize as a "vortex slush" with poly crystalline structure within the layers and interlayer 
correlations extending to about twenty layers. This is in contrast with the pinned vortex liquid 
phase into which the Bragg glass was found to melt, using the same methods, in the case of a 
large concentration of weak pinning centers: that phase was amorphous with very little interlayer 
correlation. The value of the second moment of the random pinning potential at which the Bragg 
glass melts for a fixed temperature is very different in the two systems. These results imply that 
the effects of random point pinning can not be described only in terms of the second moment of the 
pinning potential, and that some of the unresolved contradictions in the literature concerning the 
nature of the low T and high s phase in this system are likely to arise from differences in the nature 
of the pinning in different samples, or from assumptions made about the pinning potential. 
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I. INTRODUCTION 



The effects of random point pining on the properties of 
vortex matter in high-temperature superconductors have 
been extensively investigated 1 in a large number of ex- 
perimental, theoretical, and numerical studies. However, 
many questions about the phase diagram of this system 
remain controversial. The existence of a topologically or- 
dered Bragg glass (BrG) phased with long-range bond- 
orientational order and power-law decay of translational 
correlations at low temperatures and low fields in sys- 
tems with weak pinning has been experimentally estab- 
lished 4 . The BrG phase melts into a disordered vortex 
liquid (VL) as the temperature is increased. It also be- 
comes unstable as the pinning strength is increased (or 
equivalently, the magnetic field is increased for a fixed 
amount of pinning) with the temperature held fixed at 
a low value. There is no consensus at present about the 
nature of the phase into which the BrG transforms un- 
der these conditions. The theoretical proposal 5 -* 6 - of the 
existence of a vortex glass (VG) with long-range coher- 
ence of the superconducting phase and divergent barriers 
for the motion of vortex lines in the limit of zero current 
has received support from many experiments 7 -' 8 -' 9 -. How- 
ever, other experimental studies 1 ^ have questioned the 
the existence of a VG phase thermodynamically distinct 
from the high-temperature VL. Different numerical stud- 
ie a n i 12 ' 13 ] 14 have also reached different conclusions about 



the existence of the VG phase. 

The possible existence of other disordered phases, dis- 
tinct from the VL and VG phases mentioned above, has 
also been suggested in several studies. One of these 
phases is the "vortex slush" (VS) - a disordered, liquid- 
like phase with no superconducting phase coherence that 
is distinguished from the VL in that the translational cor- 
relation length in this phase is substantially larger than 
that in the VL. The existence of this phase in supercon- 
ductors with random point pinning has been predicted 
theoreticall y 15 ' 16 and several experimental studie s 17 ' 18 ' 19 
have presented evidence for the occurrence of a first-order 
transition between VL and VS phases (followed by a con- 
tinuous transition from the VS to the VG phase) as the 
temperature is decreased at fixed magnetic field. A nu- 
merical study" has also provided evidence for the exis- 
tence of the VS phase, although the validity of this ev- 
idence for macroscopic systems has been questioned in 
later work 2 ^. 

Another glassy phase proposed to exist in both high- 
temperature and conventional superconductors with ran- 
dom point pinning is the "multidomain glass" (MG) 21 ' 22 
in which the vortices are supposed to form crystalline 
domains with typical size substantially larger than the 
translational correlation length in the VL phase. In this 
MG phase, the vortex positions do not exhibit any long- 
range translational or orientational order and there is no 
superconducting phase coherence. It has been argued 
21 ' 22 that a variety of interesting "glassy" behavior ob- 
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served^ near the melting transition of the BrG phase 
may be explained by assuming that a small sliver of MG 
phase exists between the BrG and VL phases, so that the 
BrG first transforms to the MG phase as the tempera- 
ture is increased, and then the MG melts into the VL at 
a slightly higher temperature. Results of recent muon- 
spin-rotation experiment a 23 i 24 have been interpreted as 
evidence for the occurrence of vortex structures similar 
to that expected in the MG phase in the disordered glassy 
phase of a high-temperature superconductor. 

There is an additional controversy about the structure 
of the disordered vortex state at low temperatures and 
strong pinning (or high magnetic field). Early decora- 
tion experiments 2 ^ on layered high-temperature super- 
conductors showed evidence for an amorphous arrange- 
ment of vortices in the disordered phase. However, more 
recent experiments^ on a low-temperature superconduc- 
tor, NbSe2, have shown the occurrence of polycrystallinc 
disordered structures. A similar disagreement exists be- 
tween the results of different simulations: while a simu- 
lation 2 ^ of the structure of a two-dimensional vortex sys- 
tem with random point pinning shows a polycrystallinc 
arrangement of the vortices, another numerical study of 
a similar system 2 ^ has found amorphous structures with- 
out well-defined crystalline regions separated by grain 
boundaries. The reason for this disagreement remains 
to be understood. 

The question we wish to address in this work is whether 
the differences (summarized above) among the results of 
different studies of the properties of vortex matter with 
random point pinning arise, at least in part, from differ- 
ences in the details of the pinning potential. It is often 
quite difficult to obtain reliable information about the 
microscopic pinning centers present in an experimentally 
studied sample. Therefore, it is possible that samples of 
the same superconductor studied in different experiments 
have different pinning potentials. This may account for 
some of the differences in the experimental results. Simi- 
larly, the pinning potentials used in different simulations 
are often quite different from one another. It is, therefore, 
important to determine the extent to which the details 
of the pinning potential arc relevant in determining the 
properties of the system. Here we tackle this issue by 
studying the thermodynamic and structural properties 
of the mixed phase of a highly anisotropic layered high- 
temperature superconductor in the presence of a small 
concentration of randomly located strong point pinning 
centers. The pinning centers on different layers are as- 
sumed to be completely uncorrelated. The primary moti- 
vation for this study is to examine the dependence of the 
phase diagram of this system and the structure of the 
disordered phase on the details of the random pinning 
potential by comparing the results obtained with previ- 
ous results 2 ^ obtained from using the same methods for 
the same system, but with a much larger concentration 
of weaker point defects. 

In theoretical studies, the random pinning potential 
v(r) is usually assumed to be a random Gaussian variable 



with zero mean, whose statistics is completely specified 
by the second moment of its distribution, (v(r)v(r')) = 
K 2 f(\r — r'|), where f(x) is a short-ranged function (its 
range is usually assumed to be of the order of the su- 
perconducting coherence length £) normalized to unity 
at x = 0, and the parameter K measures the strength 
of pinning. For pinning due to a collection of randomly 
distributed point pinning centers, K is proportional to 
the product of the depth s of the pinning potential of 
an individual pinning center and the square root of the 
concentration c of pinning centers. Thus, two systems, 
one with small s and large c (i.e. a large concentration of 
weak pinning centers) and a second one with large s, but 
small c (i.e. a small concentration of strong pinning cen- 
ters) may have very similar values of the parameter K, so 
that an analytical treatment based on the assumption of 
Gaussian randomness would predict the same behavior 
for the two systems. Physically, however, it is not al all 
obvious that the properties of these two systems would 
be very similar. The main objective of our study is to 
examine this question. 

Our study proceeds by numerical minimization of a 
discretized version of the Ramakrishnan-Yussouff (RY) 30 
free-energy functional for the system of pancake vortices 
in a highly anisotropic layered high-tcmperature super- 
conductor, in a magnetic field normal to the layers. Dif- 
ferent phases of the system correspond to different local 
minima of the free energy in this mean-field description, 
and a crossing of the free energies of two different min- 
ima represents a first-order phase transition. Information 
about the structure of different phases is obtained from 
various correlation functions of the density distributions 
at the corresponding free-energy minima. We have previ- 
ously used this method to study the thermodynamic and 
structural properties of this system in the presence of a 
small concentration of strong columnar pinning centers 
perpendicular to the layer o 31 i 32 ' 33 , and also a large con- 
centration of randomly placed weak point pinning cen- 
ters 2 ^. We found the phase diagrams in the two cases to 
be quite different. In the system with a small concentra- 
tion of strong columnar pins, the low-temperature BrG 
phase exhibited a two-step melting, via an intermediate 
polycrystalline Bose glass phase, to the VL. This is quali- 
tatively similar to the predictio n 21 ! 22 of two-step melting 
(via an intermediate MG phase) of the BrG phase in sys- 
tems with random point pinning. In contrast, we found a 
single first-order melting of the BrG phase in the system 
with a large concentration of weak random point pin- 
ning centers. The structure of the disordered phase was 
amorphous in this case. 

The pinning potentials in these two cases differ in two 
important ways. First, the pinning potential on differ- 
ent layers are perfectly correlated for columnar pins and 
completely uncorrelated for point pinning. Second, the 
system with columnar pins has a small concentration of 
strong pining centers, whereas the one with point pin- 
ning has a large concentration of weak pinning centers. 
It is not clear which one of these two differences is the pri- 
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mary reason for the observed differences in the properties 
of these two systems. More important is the possibility 
that, in the case of point pinning alone, differences in the 
concentration and strength of the pinning centers can ac- 
count for the apparently conflicting behaviors discussed 
above. To shed light on these issues, we consider here a 
system with a small concentration of strong pinning cen- 
ters whose positions on different layers are completely 
uncorrelated. The values of the areal concentration c of 
pinning centers on each layer and the strength s of the 
pinning potential of each pinning center are chosen so 
as to includ e th e region where the values of the param- 
eter K oc V s 2 c are similar to those for the system with 
large c and small s studied earlier—. Thus, a compar- 
ison of the results obtained here with those of Ref. 
can provide useful information about the dependence of 
the properties of disordered vortex matter on the details 
of the random pinning potential. The model of pinning 
considered here may be appropriate for films of high- 
temperature superconductors in which meandering lines 
of dislocation o 34 ' 35 or artificially introduced material de- 
fects 5<> act as strong columnar pinning centers that are 
not perfectly correlated across layers. 

We consider a fixed value of the magnetic induction B 
— 2kG and study the phase behavior for different values 
of the strength s of the pinning potential, keeping the 
low pin concentration c fixed (a discussion of how the 
behavior is expected to depend on the value of B may be 
found in Ref. |29I ). For relatively small values of s, we find 
a topologically ordered BrG phase at low temperatures, 
which melts into a VL via a first-order transition as T is 
increased. As s is increased, the temperature at which 
this melting transition occurs decreases, and the transi- 
tion line in the (T — s) plane tends to become parallel 
to the T-axis for large s, indicating that the BrG phase 
does not exist if the pin strength exceeds a critical value. 
We find that the disordered phase into which the BrG 
transforms as s is increased at low T is glass-like: the 
vortices are strongly localized at points that do not form 
a structure with long-range translational or orientational 
correlation. However, this phase transforms continuously 
into the high-temperature VL as T is increased at con- 
stant s: we do not find any evidence for a phase transition 
between this glassy state and the VL. 

Although the general features of the phase diagram are 
qualitatively similar to those found in our earlier study, 29 
several other features found here are substantially differ- 
ent. The structure of the glassy state at large s and low T 
is quite different from that of the disordered phase found 
in Ref. [29i In the glassy phase of the present system, 
the vortices on each layer form a polycrystalline struc- 
ture with crystalline domains separated by grain bound- 
aries, in contrast to the amorphous structure found in 
Ref. [29|. Also, the vortex positions on different layers are 
now significantly correlated, with a correlation length of 
the order of 20 layer spacings. This is very different from 



the previously founo— glassy state, which does not ex- 
hibit any interlayer correlation in the vortex positions. 
The glassy state found here has the characteristics of the 
VS and MG phases mentioned above. For this reason, we 
classify this state as a VS, whereas the glassy state found 
in Ref. |29| was identified as a pinned vortex liquid. An- 
other difference in the properties of these two systems lies 
in the location of a crossover line that separates the re- 
gion in the (T—s) plane in which the VS is "glassy" in the 
sense that the peak vortex densities are high from that in 
which these densities are more liquid-like. This crossover 
line is obtained using a criterion^ 3 * based on percolation 
of liquid-like regions. The location of this crossover line 
in relation to the first order line in the (T — s) plane 
bounding the BrG phase is quite different in the two sys- 
tems. We also find that the rms value of the random 
pinning potential at which the BrG melts in the present 
system is about an order of magnitude smaller than that 
at which the BrG to pinned vortex liquid transition oc- 
curs in the system of Ref. at the same temperature. 
These differences illustrate the importance of the details 
of the pinning potential in determining the properties 
of disordered vortex matter, and suggest that the pre- 
dictions of analytic calculations in which the effects of 
random pinning are assumed to be described completely 
by the second moment of the distribution of the pinning 
potential may not be quantitatively correct. 

The rest of this paper is organized as follows. The 
model considered here and the method of calculation are 
described in section |TTJ In section IIII1 we describe in de- 
tail the results of our calculations and compare these re- 
sults with those obtained in 

Refill. 

Section IIVI contains 
a summary of our main results and a few concluding re- 
marks. 



II. MODEL AND METHODS 

We consider the system of pancake vortices in a 
highly anisotropic layered superconductor with vanish- 
ingly small Josephson coupling between layers in a mag- 
netic field perpendicular to the layers. In this system, 
the pancake vortices on different layers are coupled only 
via their electromagnetic interaction. Our starting point 
is the expression for the free energy of the vortex system 
as a functional of the time-averaged areal density of pan- 
cake vortices /9 ra (r). Here r is a two dimensional vector 
and the discrete index n numbers the layers. We write 

F[p] = F RY [p] + F p [p] (2.1) 

where the first term in the right-hand side is the free en- 
ergy in the absence of pinning, while the second includes 
the pinning effects. For the first term we use the RY 
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0F RY [p] =J2 J d 2 r{ Pn (r)ln(p n (r)/po) - Sp n (r)}(l/2)J2Y, J £r J d 2 r' C mn (\r - r'\)S Pm (r)5p n (r'), (2.2) 



r 



where (3 is the inverse temperature and Sp n (r) = p n {^) — 
Po the deviation of p n (r) from its average value po, the 
density of the uniform liquid (po = B/^>q where B is 
the magnetic induction and <&o the superconducting flux 
quantum). We have taken our zero of the free energy 
at its uniform liquid value. C mrl (r) denotes the direct 
pair correlation function of the uniform vortex liquid at 
density po- This static correlation function is assumed 
known and contains the required information about the 
interactions. We have taken for C mn {r) the expression 
obtained from the hypernetted chain approximation in 
Ref.l37l We use parameter values appropriate to BSCCO, 
taking the same numerical values as in Ref. [29l . 

For the second (pinning) term in Eq. (|2.1j) we write: 



^W=E / drVP(r)[ Pn (r)-p ], 



(2.3) 



where V^(r) is the pinning potential at point r on layer 
n. The pinning potential is assumed to be produced by 
random atomic scale point defects. We characterize the 
concentration of these point defects by the fraction c of 
the cells in the underlying BSCCO crystal lattice (which 
we take to be orthorhombic with in-plane lattice spacing 
do = AA and interlayer distance d = 15A) that contain 
defects. 

To study the properties of the vortex system, we dis- 
cretize the free energy by introducing variables p n ,i in 
a triangular computational lattice of size TV 2 x Nl with 
periodic boundary conditions, where Nl is the number 
of layers. Here p n ,i = Aop n (ri) where Aq is the area of 
the in-plane unit cell of the computational lattice and 
Yi denotes the location of the ith computational lattice 
point in the nth layer. We then minimize the free energy 
functional with respect to the variables p n> i using pre- 
viously described procedures ! 38 ' 39 All of the results pre- 
sented here are for computational lattices with N = 256 
and N L = 128. We work at constant field B = 0.2T. We 
take the computational lattice constant to be h = a/16, 
where a denotes the equilibrium value 38 of the lattice con- 
stant of the vortex lattice of the pure system for B = 0.2T 
near its melting transition. With these choices, the num- 
ber of pancake vortices per layer is N v = 256. All lengths 
are measured in units of ao defined by irpoaQ = 1. 

To determine the pinning potential V Ul i associated 
with each computational cell, we generate in each layer 
a random set of Nd points, where Nd is the number 
of defects per layer, Nd — cN 2 Aq/6Jq. We then have 
V nt i = s(m nt i — Nd/N 2 ) where s is the pinning poten- 
tial strength, which we will give in degrees K, and m n ^ 
the number of defect points associated with the (n, i) 
computational site. We consider here relatively large 



values of s/T and small values of c (see below). In 
that case the most frequent value of the integer m, n ^ 
is zero and it rarely exceeds unity. With these conven- 
tions the average of V n ^ is zero while its fluctuations are 
(Vn,iV n 'j) — 268.3s 2 cd nin 'Sij. In most of our calcula- 
tions, the concentration c of pinning defects was fixed 
at a value that corresponds to having Nd = 24 pins per 
layer, but we have also checked that results for 32 and 
64 pins per layer are qualitatively the same and quan- 
titatively very similar. The pin concentration of 24 per 
layer corresponds, in physical units, to having 1.36 x 10 -4 
pinning defects for every hundred unit cells in the under- 
lying crystal lattice. This is in contrast to the value used 
in Ref. |29| which was one defect per hundred unit cells 
of the crystal lattice. Thus the value of c used here is 
about four orders of magnitude smaller. In contrast we 
use values of the strength s one or two orders of mag- 
nitude larger, so that the product s 2 c that determines 
the second moment of the random pinning potential is 
comparable in the two cases. 



III. RESULTS 

In this section we discuss the results obtained using 
the methods described in the preceding section. Tests 
of the accuracy of our numerical procedures were exten- 
sively reported in Ref. [29| and need not be reported again 
here. Different phases, corresponding to different vor- 
tex density structures, are obtained by starting the min- 
imization process with different initial conditions. For a 
new pin configuration, one can start either with uniform 
conditions (all variables p n ^ set to their average value) 
or with "crystal like" initial conditions where the initial 
values of the p n .i variables are set by minimizing— the 
pinning energy of the equilibrium crystal configuration in 
the absence of pinning. The first set of initial conditions 
leads to a disordered configuration provided that such a 
configuration exists, at the values of s and T considered, 
as a local free energy minimum, that is, either a stable 
or a metastable configuration. Similarly, the second set 
of initial conditions leads to an ordered state of the BrG 
type, if such a state is at least metastable. Once a local 
minimum configuration has been obtained at a certain s 
and T, it can be used as the initial condition for a run 
at nearby values of these parameters. This procedure is 
efficient since convergence is then fast unless in chang- 
ing the parameter values the boundary of stability of the 
phase in question is crossed. We have used this method 
to carry out various annealing and thermal cycling proce- 
dures 3 ^ 3 - that ensure that the free-energy minima con- 
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FIG. 1: (Color online) Peak positions at s = 120A', T = 
13.4K for 20 consecutive layers in the disordered phase (top 
panel) and at s = 120K, T = 15.2K for all 128 layers in the 
ordered phase (lower panel). 



FIG. 2: (Color online) S(k) at T = 17.2if and s = 60K . Top 
panel: disordered phase. Bottom panel: ordered phase. 



sidered in our studies represent low-lying local minima 
of the free energy. Results for different random pin con- 
figurations at the same values of c and s are extremely 
consistent and we have found it sufficient to average over 
three different configurations in order to obtain quanti- 
ties such as the location of phase boundaries. 



A. Structure of free-energy minima 

The structure of the free energy minima can be deter- 
mined by studying different correlation functions, that 
are readily obtained from the p n .i at that local minimum, 
and by direct visualization, fn Fig Q] we consider the op- 
tion of directly visualizing the vortex lattice. From the 
values of p n ,i, we extract the average vortex positions by 
locating a vortex at site i of layer n if the value of p n ^ at 
such a site exceeds all values of p n j for any site j within 



radius a/2 oi site i. In both panels of Fig[TJ a single dot is 
plotted corresponding to the position of a vortex at any 
layer. In the top panel which is for the disordered phase 
at s — 120K , T = 13AK, this is done for twenty con- 
secutive layers, chosen at random, while in the bottom 
panel, which is for the ordered phase at s = 120K and 
T = 15.2K, this is done for all 128 layers. The difference 
between the two phases can be seen clearly from these 
plots. The ordered phase exhibits a nearly crystalline ar- 
rangement of the vortices, while no such order is visible 
in the plot for the disordered phase. If all 128 layers were 
included in the top panel, the plot would be nearly filled 
by the dots. Thus, while there are, as we will see, some 
interlayer correlations in the vortex positions in the dis- 
ordered phase, their range is much shorter than that in 
the ordered one. 

Turning now to the correlation functions, we show in 
Fig. [3 the two-dimensional static structure factor S(k) 
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FIG. 3: (Color online) The real space correlations g(r,n) as 
defined in the text, at s = WOK, T = 16.0K. The lines simply 
connect the data points. Top panel: results for the disordered 
phase for n = 0, 1,9 (from top to bottom at r — 0). Bottom 
panel: ordered phase, n — ((red) plus signs), n = 9 ((green) 
crosses) . 



defined as: 

|p(k,fc z = 0)| 2 
* (k) = N V N L ' (3 ' 1} 

where ,o(k, fc z ) is the Fourier transform of with k 
being the two-dimensional wave- vector in the layer plane 
and k z the wave-vector in the direction normal to the 
layers. With the normalization chosen (the total num- 
ber of pancake vortices) S should be of order unity for 
a disordered state. Its maximum possible value is N v Nl 
itself, which equals 32768 in our samples. Both panels 
of Fig H correspond to s = 60K, T = 17. 2K where both 
phases are locally stable. The top panel corresponds to 
the disordered phase and the lower panel to the ordered 
phase which is the globally stable one, as we will see be- 
low, at these values of (T,s). Noting the very different 
vertical scales and the well-defined hexagonal pattern in 
the lower panel, a clear difference between the two phases 
in the degree of order in the transverse planes becomes 
quite obvious. This is confirmed by the results for the 
correlations in real space, shown in Fig [31 which pertain 



also to the out-of-plane correlations. We define g(r,n) 
as the angularly averaged Fourier transform of <S(k, k z ), 
normalized by p\. Here r is the in-plane distance and n 
indexes the distance between planes. Thus, g(r, 0) is sim- 
ply the in-plane, angularly averaged correlation function 
g(r). We plot, in Fig 02 g(r, n) at s = 100K, T = 16.0K. 
At these values the ordered phase is globally stable and 
the disordered one, metastable. In the top panel (dis- 
ordered phase) we plot the cases n — 0,1,9. We see 
that the decay of g(r, n) as a function of either r or n is 
rather fast, particularly the latter. On the bottom panel, 
for the ordered phase, where we plot only the two values 
n = 0, 9, the situation is very different: the decay with 
r is markedly slower but the most remarkable feature is 
that the n = and n — 9 results are nearly identical, in- 
dicating long-range correlations in the direction normal 
to the planes. 



We now return to the vortex lattice, as in Fig. [TJ 
and re-examine it via a Voronoi construction. We re- 
call that the Voronoi construction is nothing but the 
usual Wigner-Seitz procedure on a perfect lattice, but 
performed in an arbitrary network. The number of sides 
of a Voronoi cell associated with a vertex of the network 
(location of a vortex in our case) represents the number 
of neighbors of that vertex. These numbers of neighbors 
are represent by a symbol (and color) code in Figure |U 
There we display Voronoi plots for one layer, chosen at 
random, of the systems whose vortex configurations are 
shown in Fig. [TJ One can see that in the ordered phase 
(bottom panel) nearly all the vortices have six neigh- 
bors and that the most common deviation is a pair of 
adjacent sites with five and seven neighbors, which cor- 
responds to a dislocation. These dislocations invariably 
occur near pinning sites (in fact, all the visible deviations 
from a perfect crystalline arrangement of the vortices are 
found near pinning sites which are almost always occu- 
pied by vortices, except in very rare occurrences of two 
pinning sites very close to each other), and they form 
tightly bound clusters with net Burgers vector equal to 
zero. The top plot is very different: there we can see 
nearly ordered domains in which the vortices have the 
regular number of six neighbors. These domains are small 
and separated by regions that are much more disordered, 
consisting mainly of sites with five and seven neighbors. 
These regions might be called grain boundaries if they 
were uniformly thinner. 



^From these plots and similar data obtained in the 
ranges 60K < s < 17QK, 13AK < T < 19.0K, we con- 
clude that we have here two different phases, separated 
by a first-order transition. The evidence given up to this 
point suffices to show that the ordered phase is a Bragg 
glass. The nature of the disordered phase requires more 
detailed characterization, which is addressed in the fol- 
lowing subsection. 
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FIG. 4: (Color online) Voronoi plots. The (black) dots denote 
six-neighbor sites, the (red) circles denote sites with seven 
neighbors, (blue) squares sites with five neighbors, and (ma- 
genta) triangles sites with eight neighbors. The (green) in- 
verted triangle in the bottom panel represents a site with 
four neighbors. The (black) plus signs locate pinning sites. 
Top and bottom panels are for the same configurations as in 

Fig.ru 



B. Vortex slush and pinned vortex liquid 

In this section, we present the results of a detailed 
characterization of the structure of the disordered phase. 
We compare and contrast these characterization results 
with those obtained for the disordered phase in our ear- 
lier study 29 of the structure of the disordered phase of 
the same vortex system in the presence of a much larger 
(1% of the underlying crystal lattice) concentration c 
of randomly placed weak (small s) pinning centers. To 
make this comparison meaningful, we particularly con- 
sider cases where despite the differences i n th e values of 
s and c, the values of the parameter /3V s 2 c in the two 
cases are comparable. As explained above, if the effects 
of random pinning were determined entirely by the value 



FIG. 5: (Color online) Peak positions of twenty randomly cho- 
sen consecutive layers at T = 18. 6K and a large concentration 
(1% of the underl ying crystal lattice) of weak (s — 4K) point 
defects, as in Ref. l29l . Compare with Fig. [1] top panel. 



of s 2 c, then the structures of the disordered phases in 
the two systems would have been very similar. As we will 
seen in detail, we find very significant differences between 
the structures of the two disordered phases: indeed, we 
show that while the disordered phase in the system with 
a large concentration of weak pinning centers has very lit- 
tle short-range order (this led us to conclude^ that this 
phase should be identified as a pinned vortex liquid), the 
disordered phase in the present system with a small con- 
centration of strong pinning centers exhibits substantial 
short-range order, both in the layer plane and in the di- 
rection normal to the layers, leading us to identify this 
phase as VS or MG. 

We begin with a qualitative comparison of the peak- 
position plots (see Fig.[l] top panel), and of the Voronoi 
plots (see Fig 21 top panel) . In Fig. [5J we have shown 
the positions of the local density peaks on randomly se- 
lected 20 consecutive layers of a disordered free-energy 
minimum obtained at T = 18. 6K for a sample 29 with the 
larger value of c and weak (s = 4K) pinning centers. The 
dots that represent the peak positions fill the sample area 
quite randomly, in contrast to the plot shown in the top 
panel of Fig. [TJ which shows a definite short-range struc- 
ture, indicating the presence of correlations among the 
positions of the local density peaks. Similarly, a typical 
Voronoi plot for a randomly chosen layer of the minimum 
for which results are given in Fig. [5] is shown in Fig. [5] 
Again, in comparison with the analogous plot shown in 
the top panel of Fig. 2J this Voronoi plot shows a much 
higher degree of disorder. In particular, the total number 
of "vortex lattice" sites with six nearest neighbors is sub- 
stantially smaller in Fig. [5l and there is no evidence for 
the presence of any crystalline domain (for this reason, 
we classified^ the structure of these minima as amor- 
phous, not polycrystalline, in our earlier study). These 
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FIG. 6: (Color online) Voronoi plot at T — 18.6K and a large 
concentration (1% of the underlying crystal lattice) of weak 
(s = 4K) point defects, as in Ref. 123 . The meaning of the 
symbols is the same as in the top panel in Fig [4] with which 
the current plot should be compared. 



FIG. 7: (Color online) The in-plane vortex correlation func- 
tion C(r) (see text) as a function of r, for the disordered 
minimum obtained in the present calculation at s = 170A', 
T = 14.6K ((red) solid line), compared with the result for 
weak point defects shown in Figure |6] ((blue) dotted line). 



two plots, therefore, suggest that the disordered minima 
obtained in the present study exhibit a much higher de- 
gree of short-range order than those found earlier for a 
system with a large concentration of weak pinning cen- 
ters. 

To make these observations more quantitative, we have 
calculated two correlation functions, C(r) and C z (n), 
that represent, respectively, the degree of correlation 
among the vortex positions (i.e. the positions of the local 
peaks of the density field) in the layer plane and across 
layers. The function C(r) is the usual angularly aver- 
aged pair-distribution function of vortex positions on the 
same layer, averaged over all the layers. It should not be 
confused with the correlation function g(r) = g(r, n = 0) 
considered in section llll Al which represents the two-point 
correlation of the time-averaged local density: therefore, 
g(r) is sensitive to both the positions and the heights 
of the local density peaks, whereas C(r) provides more 
accurate information about the correlation in the peak 
positions because it is insensitive to the heights of the 
local density peaks. The function C z (n) is defined as fol- 
lows: For each peak position on a layer, we calculate the 
number of peaks on another layer, separated by n layer 
spacings from the original one, that lie within a small in- 
plane distance ro from it. This number is then averaged 
over all the peak positions on all the layers and divided 
by tttqPo, the average number of vortices in a circular 
area of radius ro , to obtain C z (n) . The values of C z (n) 
are found to be insensitive to the choice of ro as long 
as ro is small compared to the average distance between 
nearest-neighbor vortices (~ 1.9a , see Fig. [3l top panel) 
on the same layer. The results reported here were ob- 
tained using ro = 0.22ao. It is clear from the definition 
of C' z (n) that this function measures the degree of align- 



ment between vortices on two layers separated by n layer 
spacings. If the vortex positions on different layers are 
completely uncorrelated, then C z (n) should be equal to 
unity for all n. Values of C z (n) > 1 indicate some degree 
of alignment between the vortex positions on different 
layers. 

Typical results for these two correlation functions are 
shown in Figs. [7] and The results for a system with 
a small concentration of strong pinning centers (present 
study) are shown for s — 170K and T = 14. 6K, whereas 
the results for a system with a large concentration of 
weak pinning centers (Ref. l29t ) are for the same parame- 
ters as in Figs. [5] and [6] These were chosen so as to make 
the values of /3V s 2 c for the two systems very similar. 
Also, the functions C(r) and C z (n) that describe corre- 
lations among average vortex positions are meaningful 
only when the local density peaks are sufficiently high 
(the vortices are strongly localized), so that the average 
vortex positions can be defined without ambiguity This 
criterion is satisfied for the parameter values for which 
results are shown in Figs. [7] and [3 

The plots in Fig. [7] clearly show that the disordered 
phase in the present (low c) system exhibits a substantial 
degree of in-plane (short-range) positional order, while 
the disordered phase in the system with a large concen- 
tration of weak pinning centers has almost no in-plane 
correlation in the vortex positions. The function C(r) 
for this low c and high s system exhibits six clear peaks 
of decreasing height, indicating the existence of in-plane 
positional correlations up to a distance of about 10ao- 
This length scale is consistent with the typical size (20-30 
vortices, see Fig. 21 top panel) of in-plane crystalline do- 
mains at the disordered free-energy minima found in the 
present work. In contrast, the C(r) for the system with 
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FIG. 8: (Color online) The out-of-plane vortex correlation 
function C' z (n) (see text) shown for the same two cases as in 
the preceding figure. The (red) circles are the high s case of 
the present work. The solid line is an exponential fit. 



a large concentration of weak pinning centers exhibits no 
prominent peaks, indicating that the disordered phase of 
this system has very little in-plane positional order and 
crystalline domains can not be meaningfully defined in 
this case. 

A similar difference in the degree of out-of-plane cor- 
relations in the vortex positions is also found in Fig. [8j 
where the function C z (n) is plotted for the two free- 
energy minima of Fig. [7] The data for the system with a 
large concentration of weak pinning centers are all close 
to 1.0, even for values of n near 1. This implies that 
the vortex positions on different layers are essentially un- 
corrected in this case. On the other hand, C z (n) for 
the present system starts at a value close to 2.0 for n 
near 1, and approaches unity only as n approaches about 
Nl/3, (the largest value of n for which C z (n) can be 
defined in our system of Nl layers with periodic bound- 
ary conditions is Nl/2). These results indicate that the 
vortices on different layers are aligned to some extent 
in the disordered phase of the present system. We find 
that the dependence of C z (n) on the layer separation 
n can be represented quite well by the functional form 
C z (n) — 1.0 + C z q exp(— n/l z ) where the "correlation 
length" l z provides a measure of the degree of out-of- 
plane alignment of the vortices. A fit of the data to this 
functional form (with C z q — 1.0, l z — 17.2) is shown by 
the solid line in Fig. [51 Similar values of l z are found for 
other values of the parameters in the disordered phase of 
the present system, indicating that the vortex positions 
on different layers are substantially correlated in the sys- 
tem with a small concentration of strong pinning centers. 

It is clear from all the results described above that 
the disordered free-energy minima found in the present 
study exhibit a substantial degree of short-range posi- 
tional order. The in-plane structure can be described 





n nno 




-U.UU4 - 




-0.006 - 




-0.008 - 


'"Hi 

> 


-u.in - 


z 




LL 


-U.U I d. - 








-U.U14 - 




-0.016 - 




-0.018 - 




-0.02 - 




18 



18.8 



FIG. 9: (Color online) Free energy crossings at s = 60K as 
a function of T. The quantity plotted is the free energy per 
vortex in units of feT. The (red) plus signs are results for the 
VS phase and the (green) x signs for the BrG phase. Straight 
lines connect the data points. 



as polycrystalline, with 20-30 vortices in each crystalline 
domain. The vortex positions on different layers are also 
correlated, with a correlation length of 15-20 layer spac- 
ings. These features are similar to those expected for the 
VS and MG phases that have been predicted to exist in 
vortex systems with random point pinning. We, there- 
fore, identify the disordered phase in the present system 
as VS. This is qualitatively different from the disordered 
phase found in our earlier study2£ of a vortex system 
with a large concentration of weak pinning centers. In 
that case, the disordered phase, which does not exhibit 
any appreciable in-plane or out-of-plane correlation in 
the vortex positions, was identified as pinned vortex liq- 
uid which should be distinguished from the VS found in 
the present work. 



C. Phase Diagram 

To determine the phase diagram one simply looks for 
free energy crossings as a function of s and T . This is 
quite straightforward. An example is shown in Figure [9l 
where we consider the free energy per vortex, in units of 
/csT, for the BrG and VS phases at constant s as a func- 
tion of temperature. In the range plotted both phases 
are locally stable. The BrG minimum becomes unsta- 
ble just beyond the right edge of the plot. The symbols 
correspond to the results for the free energy and they 
are connected by straight lines. Error bars from statisti- 
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cal averaging over different pin configurations are smaller 
than the symbols. The first order transition between the 
two phases is easily located from such data. A similar 
first order transition from the BrG to the VS phase is 
found as the pinning strength s is increased while keeping 
T constant at a low value. The BrG minimum becomes 
unstable as s or T is increased slightly above the value 
at which the transition to the VS phase occurs. The VS 
minimum, on the other hand, remains locally stable for 
all the values of s and T considered here. 

For relatively small values of s and large T, the val- 
ues of the density at the local peaks of the density field 
in the VS phase do not greatly exceed (see below) the 
average liquid density value po except for peaks at the 
pinning centers. In that sense, the VS phase might be 
said to be liquid-like. On the other hand, for large s and 
low T, these local peak density values are substantially 
larger than po, indicating strongly localized vortices. A 
VS minimum with these properties should be viewed as 
solid- like (glass). However, for the parameter range ex- 
plored here, we do not find any evidence for a transition 
within the VS state between the high T-low s and the 
low T-high s structures. A liquid-like minimum obtained 
by starting the free-energy minimization from disordered 
initial conditions at a high T evolves continuously to a 
glassy VS minimum as it is "followed" to a low tempera- 
ture by reducing T in small steps and performing the free- 
energy minimization at each new temperature with the 
minimum obtained at the previous temperature as the 
initial state. Similarly, a glassy VS minimum obtained 
at low T and large s by starting the minimization from 
disordered initial conditions (or by starting from a BrG 
minimum and increasing s to a value at which the BrG 
minimum becomes unstable and a disordered minimum 
is found) evolves continuously to the high-temperature 
liquid-like structure as T is increased in small steps. 

Thus, our results do not show the first-order VL to 
VS transition found in some experiment s 17 ' 18 ' 19 and in 
a simulation 14 . However, these experiments and simu- 
lation show that the line of first-order VL to VS tran- 
sitions ends at a critical point as the pinning strength 
(in the simulation) or the magnetic field (in the experi- 
ments) is increased. It is likely that the parameters used 
in our study correspond to values beyond this critical 
point, so that no VL to VS transition is found. In a 
stud y 31 ' 32 ' 33 of the same vortex system in the presence 
of a small concentration of strong columnar pins (with 
values of s comparable to those used here) perpendicu- 
lar to the layers, we found a first-order transition from 
the high-temperature VL to a polycrystalline Bose glass 
phase as T was decreased, but only when the ratio of the 
number of pinning centers to the number of vortices per 
layer was less than about 1/32. The value of this ratio in 
the present study (24/256 = 3/32) is larger than 1/32. 
This suggests that a first-order VL to VS transition may 
occur in the present system for a smaller concentration of 
pinning centers. It is, however, difficult to test this possi- 
bility numerically. The typical size of crystalline domains 




FIG. 10: (Color online) Percolation probability plotted as a 
function of the fraction / of liquid like sites. The (red) dots 
are the data, as explained in the text. The (blue) dashed line 
is along the value of 1/2. 



in a VS state will increase as the concentration of pin- 
ning centers is reduced. A multi-domain structure can be 
distinguished from a crystalline one in a numerical study 
only if the sample size is substantially larger than the typ- 
ical domain size. So, to study polycrystalline free-energy 
minima in the present system for a much smaller con- 
centration of pinning centers, one would require samples 
containing a few thousand vortices per layer. Samples of 
such size could be studied for columnar pins because the 
problem is then effectively two-dimensional, but a similar 
study for the present case (where the problem is three- 
dimensional) would be computationally intractable. 

A crossover line between liquid-like and glassy behav- 
iors of the VS phase can be drawn in the (T — s ) pl ane 
by using a criterion based on percolation. In Ref. l33l . we 
found that a local peak of the density field may be classi- 
fied as "liquid-like" (the corresponding vortex is weakly 
localized) if p pea k , the value of the density variable at the 
peak, does not exceed 3po- Otherwise it is characterized 
as "solid-like" , representing a strongly localized vortex. 
One can then readily examine whether or not, at a certain 
s and T, the regions with Ppeak 

< 3po percolate across 
any of the layers in a sample. If such regions percolate in 
a majority of the layers, then the corresponding VS min- 
imum may be classified as "liquid-like"; otherwise, the 
minimum should be called "glassy" . Similar percolation 
criteria have been used in other studie s 40 ' 41 to differen- 
tiate between liquid and glassy phases. We have done 
this for all values of s and T studied. The probability 
that liquid-like sites percolate across a layer depends on 
the fraction / of such sites. This dependence is shown 
in Fig. [10l where we see that the percolation probability 
is 1/2 very close to / = 1/2. We have therefore adopted 
the criterion that percolation of liquid-like regions occurs 
at / = 1/2, which is also the percolation threshold for 
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random site percolation on a triangular lattice. 

The phase diagram of the system in the (T — s) plane 
is shown in Fig. [ITJ top panel. The plus signs connected 
by solid red segments locate the line of first order transi- 
tions separating, as explained above, the BrG and the 
VS phases. We have also indicated in the panel the 
percolation line for high s and low T - percolation of 
liquid-like regions in the VS phase occurs as this line is 
crossed from the left or from above. A few features of 
the phase diagram are worth pointing out. The temper- 
ature at which the melting transition of the BrG phase 
occurs is quite insensitive to the value of s for relatively 
small values of s. As s is increased further, the transition 
temperature decreases and the phase boundary tends to 
become parallel to the temperature axis for large s. The 
decrease in the transition temperature with increasing 
s is in qualitative agreement with the prediction of ana- 
lytic calculations i 42 ! 43 ' 44 The BrG phase exists only if the 
pinning strength is sufficiently small. The shape of the 
BrG-VS phase boundary in our phase diagram is similar 
to that found in the simulation of Ref. [lj in which the 
ratio of the number of pinning centers to the number of 
vortex lines is similar to the value used here. The loca- 
tion of the percolation line relative to the BrG-VS phase 
boundary is also similar to that of the VS-VG transition 
line found in Ref. [Til . This, however, may be a coinci- 
dence, since the percolation line in our phase diagram 
does not correspond to a true phase transition, whereas 
the VS-VG transition line in Ref. [13 represents a phase 
transition signaled by the appearance of superconducting 
phase coherence which can not be studied in our calcu- 
lation. 

For comparison we have included in Fig. [TTJ (lower 
panel) the phase diagram obtained 29 for a high concen- 
tration of weak pinning centers, including also the per- 
colation line, which was not reported earlier. The gen- 
eral shapes of the phase boundaries obtained in the two 
cases are similar. However, a quantitative comparison 
shows that a low concentration of strong pinning cen- 
ters considered in the present study is much more ef- 
fective in destroying the nearly crystalline order of the 
BrG phase than a large concentration of weak pinning 
centers considered previously. For example, in the phase 
diagram shown in the top panel of Fig. 111! the transi- 
tion of the BrG phase to the VS phase at T = 17 K oc- 
c urs at s ~ 100A, for which the rms pinning potential is 
V268.3s 2 c = 1.9K (see Section [TTJ) . In contrast, the tran- 
sition at T = 17K in the phase diagram in the bottom 
panel of Fig. [TT] occurs near s = 11 A, for which the rms 
value of the pinning potential is 18. OA". Thus, the rms 
value of the random pinning potential at which the BrG 
phase undergoes a transition to a disordered phase at a 
fixed temperature in a system with a large concentration 
of weak pinning centers is about 10 times larger than that 
in a system with a small concentration of strong pinning 
centers. This implies that the kind of pinning considered 
in the present study is much "stronger" (i.e. more ef- 
fective in destroying BrG order) than that considered in 




FIG. 11: (Color online) Phase diagram (top panel) for the 
present case. The (red) plus signs are the results for the 
first order transition between the BrG and VS phases. The 
points are connected by straight line segments. The (green) 
x signs connected by a dotted line denote the points at which 
percolation of liquid-like region occurs in the VS phase. For 
comparison, the bottom panel shows the phase diagram at a 
high concentration of weak pinning centers (Ref. [5^) including 
the percolation line. 
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Ref. [H. 

The comparison between the two phase diagrams also 
points out a serious flaw in analytic calculations in which 
the random pinning potential is assumed to be Gaussian 
and its effects are assumed to be determined completely 
by its rms value (second moment). Our results show that 
the details of the pinning potential, not just its second 
moment, are very important in determining the phase 
behavior: the value of the second moment of the pinning 
potential at which the BrG phase undergoes a transition 
to a disordered phase at a fixed temperature can vary 
by orders of magnitude, depending on the details of the 
pinning potential. Thus, the predictions of analytic cal- 
culation o 42 ' 43 ' 44 ' 45 ' 46 of the phase diagram of supercon- 
ductors with random pinning, in which the effects of the 
pinning potential are assumed to be determined by its 
second moment only, can not be quantitatively accurate 
- the details of the pinning potential have to be taken 
into account for an accurate theoretical determination of 
the locations of the phase boundaries. 

Fig, mi also shows that the location of the percolation 
line relative to that of the transition line of the BrG phase 
is very different in the two cases. In the top panel, the 
percolation line intersects the BrG-VS transition line at 
a relatively low temperature where the latter is almost 
parallel to the T-axis, whereas in the bottom panel, the 
intersection is near the high-temperature, "vertical" part 
of the boundary of the B rG phase. The value of the di- 
mensionless quantity /3V s 2 c, which measures the ratio 
of the rms value of the pinning potential to the thermal 
energy is, along the percolation line, nearly the same in 
the two cases, indicating that the rms value of the pin- 
ning potential provides a fairly accurate account of the 
effectiveness of random pinning in localizing the vortices. 
The difference in the location of the percolation line rela- 
tive to the phase boundary of the BrG phase again arises 
because the rms value of the pinning potential is not, by 
itself, a good indicator of the effectiveness of the pinning 
potential in destroying BrG order. 

IV. SUMMARY AND DISCUSSION 

We have considered here a layered superconductor with 
a small concentration c of strong random point pinning 
centers and a fixed magnetic induction, normal to the lay- 
ers. We have described the phase diagram in the (T — ,s) 
plane, where s is the pinning strength. We find that the 
stable phase at sufficiently low values of T and s is clearly 
a Bragg glass (BrG), as in the high c case^ The melting 
line in the (T — s) plane at low T is nearly horizontal. 
Upon increasing s or T the BrG melts via a first order 
transition into a disordered phase. We find that this 
phase is clearly polycrystalline (see top panels in Figs [1] 
and 2]) with well defined crystalline domains separated 



by domain walls. It has characteristics of the VS phase 
discussed in the hterature ) 14 ' 15 i 16 i 17 i 18 i 19 and of the also 
propose d 21 ' 22 MG phase. We therefore characterize this 
phase (see Fig. [TTJ and the text) as a VS. There is a very 
clear difference between this and the high c, low s case 
as to the nature of this disordered phase. For large c and 
low s, this phase was found 2 ^ to be an amorphous, pinned 
VL state, with almost no correlations in the layer plane 
or between different layers (see Figs. [6] and [7]). A key 
result of this paper is this contrast between the behavior 
of the system studied here and that of the correspond- 
ing system 2 ^ with a much larger concentration of weak 
pinning centers. 

The phase diagrams of these two systems in the (T — s) 
plane are qualitatively similar, but very different quanti- 
tatively. The value of the second moment of the random 
pinning potential at which the transition of the BrG to 
the VS occurs in the present system is much smaller than 
that at the BrG to pinned VL transition in Ref. l20j at the 
same temperature. Also, the location of the percolation 
line relative to the phase boundary of the BrG phase 
is quite different in the two systems. These results im- 
ply that the notion that the parameter j3\/ s 2 c suffices 
to characterize the effects of the random point pinning 
potential (in other words, that it is enough to specify 
its second moment) is too simplistic and may be quan- 
titatively inaccurate. The assumption of Gaussian ran- 
domness, which implies that only the value of the second 
moment of the distribution of the random pinning po- 
tential is relevant, is valid only in a rather limited sense, 
namely only when the pinning potential seen by a vor- 
tex is the sum of contributions from a large number of 
weak pinning centers. The validity of the Gaussian as- 
sumption does not extend further than that, as we have 
proved here. Although the limited validity of this as- 
sumption was mentioned in a few existing studies , 42 ' 43 
which did not however quantify its degree of inadequacy, 
it has been generally ignored. 

Our results suggest that many of the experimental and 
theoretical controversies found in the literature (which we 
have reviewed in our Introduction) about the behavior of 
this system originate in differences in the insufficiently 
characterized nature of the point pinning. This con- 
clusion, besides its obvious theoretical import, has also 
strong practical implications since a good understanding 
of the effects of pinning is the key to obtaining the larger 
values of the critical current that are needed for the high 
temperature superconductors to fulfill their promise as 
potentially useful materials. 
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